% File: sc_fde_ib.m
% -----------------
% This script implements baseline ib-dfe rx

H_tmp  = squeeze(H).';
H_tmp  = reshape(H_tmp,[],nt);
Var_S  = ones(nsubf,1);
s_mean = zeros(nsymb,nsubf);
for p = 1: nsubf
    N_IB = 5; % iteration limit of ib-dfe
    n_ib = 0; % reset iteration counter
    while (n_ib<N_IB)
        Fib       = conj(H_tmp)./(abs(H_tmp).^2*Var_S(p)+noisePow_norm);
        gamma     = mean(Fib.*H_tmp,1);
        var_dist  = abs(1/gamma-Var_S(p)); % output mse
        % -----------------------------------------------------------------
        % test point - mse calculate
        % -----------------------------------------------------------------
        mse_buffer(n,n_ib+1,p) = mse_buffer(n,n_ib+1,p) + var_dist;
        % -----------------------------------------------------------------
        Fib = Fib./gamma;
        S_mean    = fft(s_mean(:,p),nfft,1);
        S_pre_det = S_mean + Fib.*(fdein(:,p)-H_tmp.*S_mean); % ib-dfe
        s_pre_det = ifft(S_pre_det,nfft,1);
        % -----------------------------------------------------------------
        % test point - normalized MSE of equalizer output
        % -----------------------------------------------------------------
        mse_test = mean(abs(s_pre_det - x_pload(:,p)).^2)/xtxPower_mean;
        mse_test_buffer(n,n_ib+1,p) = mse_test_buffer(n,n_ib+1,p) + mse_test;
        % -----------------------------------------------------------------
%         x_pre_det = repmat(s_pre_det,[1,M_mod]);
%         app       = exp(-1*(abs(x_pre_det-symb_sample).^2)/(var_dist*xtxPower_mean)); % symbol wise a.p.p.
%         app_norm  = repmat(1./sum(app,2),[1,M_mod]);
% %         % -----------------------------------------------------------------
% %         % test point - app
% %         % -----------------------------------------------------------------
% %         app_test       = exp(-1*(abs(x_pre_det-symb_sample).^2)/(2*var_dist*xtxPower_mean))/(sqrt(2*pi*var_dist*xtxPower_mean));
% %         app_norm_test  = repmat(1./sum(app_test,2),[1,M_mod]);
% % %         app      = app_test;
% % %         app_norm = app_norm_test;
% %         % -----------------------------------------------------------------
%         app       = app_norm.*app;
%         s_mean(:,p) = sum(symb_sample.*app,2); % mean of app detection
%         pow_s     = sum(abs(symb_sample).^2.*app,2); % mean square of app detection
%         var_s     = pow_s - abs(s_mean(:,p)).^2; % variance of app detection
%         Var_S(p)  = mean(var_s,1)/xtxPower_mean;
%         % -----------------------------------------------------------------
%         % test point - normalized VAR of soft estimator output
%         % -----------------------------------------------------------------
%         var_buffer(n,n_ib+1,p) = var_buffer(n,n_ib+1,p) + Var_S(p);
%         var_test = mean(abs(s_mean(:,p) - x_pload(:,p)).^2)/xtxPower_mean;
%         var_test_buffer(n,n_ib+1,p) = var_test_buffer(n,n_ib+1,p) + var_test;
%         % -----------------------------------------------------------------
        % -----------------------------------------------------------------
        % test point - soft detection I/Q
        % -----------------------------------------------------------------
        s_pre_det_real = real(s_pre_det);
        s_pre_det_real_soft_pos = exp(-1*(s_pre_det_real-1).^2/(var_dist*xtxPower_mean));
        s_pre_det_real_soft_neg = exp(-1*(s_pre_det_real+1).^2/(var_dist*xtxPower_mean));
        s_mean_real    = (s_pre_det_real_soft_pos - s_pre_det_real_soft_neg)./(s_pre_det_real_soft_pos + s_pre_det_real_soft_neg);
        s_pre_det_imag = imag(s_pre_det);
        s_pre_det_imag_soft_pos = exp(-1*(s_pre_det_imag-1).^2/(var_dist*xtxPower_mean));
        s_pre_det_imag_soft_neg = exp(-1*(s_pre_det_imag+1).^2/(var_dist*xtxPower_mean));
        s_mean_imag    = (s_pre_det_imag_soft_pos - s_pre_det_imag_soft_neg)./(s_pre_det_imag_soft_pos + s_pre_det_imag_soft_neg);
        s_mean_test    = s_mean_real + 1i*s_mean_imag;
        real_mean      = mean(s_mean_real);
        real_var       = mean(s_mean_real.^2 - real_mean.^2);
        imag_mean      = mean(s_mean_imag);
        imag_var       = mean(s_mean_imag.^2 - imag_mean.^2);
        s_mean_var     = real_var + imag_var;
        var_soft_test1 = abs(mean(2 - abs(s_mean_test).^2))/xtxPower_mean;
        %var_soft_test1 = 1 - s_mean_var/xtxPower_mean;
%         var_soft_test1 = 10^(a_cf * var_dist^b_cf + c_cf);
        s_mean(:,p)    = s_mean_test;
        Var_S(p)       = var_soft_test1;
        % -----------------------------------------------------------------
        % -----------------------------------------------------------------
        % test point - normalized VAR of soft estimator output
        % -----------------------------------------------------------------
        var_buffer(n,n_ib+1,p) = var_buffer(n,n_ib+1,p) + Var_S(p);
        var_test = mean(abs(s_mean_test - x_pload(:,p)).^2)/xtxPower_mean;
        var_test_buffer(n,n_ib+1,p) = var_test_buffer(n,n_ib+1,p) + var_test;
        % -----------------------------------------------------------------
        n_ib = n_ib + 1;
    end
end
z = demodulate(modem.qamdemod(M_mod),s_mean);

% End of script